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Abstract 

Ground-based arrays of imaging atmospheric Cherenkov telescopes have emerged as 
the most sensitive 7-ray detectors in the energy range of about 100 GeV and above. 
The strengths of these arrays are a very large effective collection area on the order 
of 10^ m'^, combined with excellent single photon angular and energy resolutions. 
The sensitivity of such detectors is limited by statistical fluctuations in the num- 
ber of Cosmic Ray initiated air showers that resemble gamma-ray air showers in 
many ways. In this paper, we study the performance of simple event reconstruction 
methods when applied to simulated data of the Very Energetic Radiation Imaging 
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Telescope Array System (VERITAS) experiment. We review methods for recon- 
structing the arrival direction and the energy of the primary photons, and examine 
means to improve on their performance. For a software threshold energy of 300 
GeV (100 GeV), the methods achieve point source angular and energy resolutions 
of crg3% = 0.1° (0.2°) and 0"g8% = 15% (22%), respectively. The main emphasis of 
the paper is the discussion of 7-hadron separation methods for the VERITAS exper- 
iment. We find that the information from several methods can be combined based 
on a likelihood ratio approach and the resulting algorithm achieves a 7-hadron sup- 
pression with a quality factor that is substantially higher than that achieved with 
the standard methods used so far. 

Key words: Gamma-Ray Observations, Data Analysis Methods, Imaging 
Atmospheric Cherenkov Telescopes 



1 Introduction 



The combination of the imaging optics and pixilated camera of the Whipple 10 
m Cherenkov telescope and the "Hillas" parameterization [1] of the air shower 
images led to the initial detection of TeV 7-rays from the Crab Nebula during 
1986 to 1988 [2]. The discovery established ground-based TeV 7-ray observa- 
tions as an exciting field at the intersection of astronomy and particle physics. 
Today, almost two decades later, TeV 7-ray astronomy made a major impact 
on our understanding of the blazar class of Active Galactic Nuclei [3,4,5] and 
is beginning to have a major impact on galactic high energy astronomy [6]. 

The hardware of ground-based Cherenkov telescopes has evolved significantly. 
Four large Cherenkov observatories have been built or are under construc- 
tion: VERITAS [7], the High Energy Stereoscopic Array (H.E.S.S.) [8], the 
Major Atmospheric Gamma Imaging Cherenkov detector (MAGIC) [9], and 
CANGAROO III [10]. All these four experiments use (or will use) several 
telescopes to detect air showers in coincidence. In the case of the VERITAS 
array, at first four, then later seven telescopes are to be located a distance of 
80 m from the central telescope. The H.E.S.S. telescope array is already fully 
operational and consists of four telescopes located at the corners of a square 
with sides of 120 m length. In CANGAROO HI, the telescopes are located 
on the corners of a parallelepiped with sides of 100 m length. The MAGIC 
experiment will use two telescopes at a distance of 80 m. The "stereoscopic" 
detection of air showers with several telescopes under different viewing angles 
suppresses background events initiated by muons passing close to individual 
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telescopes and simplifies the reconstruction of the air shower parameters (ar- 
rival direction and energy of the primary particle). While the Whipple 10 m 
telescope used a camera of 37 pixels (0.5° pixel pitch) two decades ago, the 
modern counterparts use reflectors of 10 m to 17 m diameter together with 
cameras containing 500 to 960 pixels (~0.15° pixel pitch). These finely pixi- 
lated cameras are read out with fast electronics. In the case of the VERITAS 
experiment, the signals are sampled at 500 MHz, making it possible to acquire 
the air shower signals within a few nsec, minimizing the contamination of the 
fast air shower pulses with night sky background photons. 

Stereoscopic observations of air showers with several Cherenkov telescopes re- 
quires a set of analysis methods to reconstruct the air shower parameters, 
and thus estimate the arrival direction as well as the type and energy of the 
primary particle. Several different approaches have been used, and we discuss 
here simple "geometrical methods" [11,12,13] and so-called "template fitting 
methods" [14,15]. The interested reader might consult [16] for a description of 
geometric methods specifically aimed at low-energy events. 

In the first approach, images are cleaned to suppress signals from night sky 
background photons and characterized with Hillas second moment parameters. 
A Cherenkov image of an air shower looks like a 2-D Gaussian distribution, or 
like a "filled ellipse" ; the second moment analysis gives the orientation of the 
major axis of the ellipse and the "width" and "length" parameters. The latter 
two parameters correspond to the root mean square (RMS) of the signal am- 
plitudes perpendicular to the major and minor axis of the image, respectively. 
Stereoscopic observation of an air shower with several telescopes under signif- 
icantly different viewing angles makes it possible to combine the information 
from the air shower images taken with all the telescopes to reconstruct the ori- 
entation of the air shower axis and its location relative to the telescopes with 
a simple geometric method. Once the air shower axis has been determined, 
information about the type of the primary particle (photon or hadron) and 
its energy can be inferred based on the Hillas parameters. 

In the second approach, a dataset of simulated photon-initiated air showers is 
used to derive the expected distribution of signal amplitudes in all telescopes 
of the experiment as function of a set of air shower parameters (e.g. orientation 
and location of the image, the energy of the primary particle, and the height 
of the shower maximum). An air shower event is then analyzed by determining 
for which parameter combination the expected signal amplitude distributions 
best agree with the observed ones. In other words, the air shower parameters 
are determined by finding the template image that resembles the observed 
image as closely as possible and minimizes a goodness-of-fit-measure like the 
X^-value or a log-likelihood parameter. The best fit gives a direct estimate of 
the air shower parameters, and the goodness-of-fit-measure of the best fit can 
be used to infer information about the likelihood that the event was initiated 



by a photon rather than a hadron. 

In this paper, we follow the first approach and evaluate its performance when 
applied to simulated VERITAS data. Direct comparison of the geometrical 
and template-fitting approaches with the HEGRA experiment [17] showed 
that the first approach can achieve an angular resolution, energy resolution 
and 7-hadron separation capability similar to the second [13,15] ^ . The geo- 
metrical approach has the important advantage that it minimizes the depen- 
dence on the input from Monte Carlo simulations and that it is computa- 
tionally less demanding. The reconstruction of the air shower axis does not 
require any Monte Carlo input at all. The 7-hadron separation method and 
the energy reconstruction method require minimum simulation input, and ex- 
perimental air shower data can be used to test that the simulations correctly 
describe the relevant air shower properties. For example, experimental 7-ray 
air shower data can be used to verify that the simulations correctly predict 
the lateral Cherenkov light distribution [18] as well as the width of the images 
as function of signal amplitude and shower axis distance [19] . Thus, while the 
template-fitting approach minimizes the statistical errors on the reconstructed 
air shower parameters, the geometrical method minimizes the systematic er- 
rors on the analysis results. We think that the robustness of the geometrical 
approach will guarantee its continued use. 

To limit the scope of the paper we have focused the discussion on the major- 
ity of the events detected with VERITAS, which means events close to the 
~100 GeV energy threshold of the experiment. While most of the methods 
presented in the following work well over the entire VERITAS energy range ^ , 
we have not specifically optimized the methods for giving the best results at 
the high-energy end. The reader is referred to [20] for a discussion of the astro- 
physics at >10 TeV energies, and experimental approaches to obtain excellent 
sensitivities in this energy range. 

Another limitation of this paper is that we study only vertically incident "low 



^ The CAT collaboration applied a template-fitting method to the data of the CAT 
telescope and reported that the method performed markedly better than a simple 
geometrical method [14]. The difference of the HEGRA and CAT experiences may 
have several reasons. It may be that the geometrical method performs relatively 
better on stereoscopic data (HEGRA) than on single telescope data (CAT). The 
cameras of the HEGRA and CAT experiments had pixels with angular diameters of 
0.25° and 0.12°, respectively. The template-fitting method may show its full power 
only when used with a finely pixilated camera. 

^ Figs. 3 and 6 show that VERITAS has a good collection efficiency and angular 
resolution at high energies; the "proton" curves in Fig. 10 show that the hadron 
suppression works well at high energies; finally, the "electron" curves in Fig. 10 
show that the cuts accept a high fraction of electromagnetic air showers even at 
high energies. 



zenith angle" air showers. Literature on the analysis of large zenith angle data 
can be found in [21,22,23] and in the references therein. 

The rest of this paper is organized as follows. In Sect. 2 we describe the shower 
and detector simulation codes used, and discuss the simulated dataset. Fur- 
thermore, we look at some raw parameter distributions, like the shower core, 
direction, and energy distributions of the events that triggered the telescope 
array. Subsequently, we briefly describe the methods used to reconstruct the 
air shower axis in Sect. 3 and the primary energy in Sect. 4. The information 
about the shower axis reconstruction and the energy estimator will be used 
by the 7-hadron separation methods. We address the suppression of Cosmic- 
Ray initiated air showers in Sect. 5. The HEGRA and H.E.S.S. experiments 
use a "scaled width parameter" described further below to separate photons 
from Cosmic Rays. Here, we evaluate several 7-hadron separation parameters. 
The general ideas behind all the 7-hadron separation methods that we will 
explore were discussed by Hillas and co-workers (see, e.g. [27]) based on very 
detailed studies of the physical properties of photon and Cosmic-Ray initiated 
air showers. Here we study specific implementations of the methods, and eval- 
uate their performance when applied to data from the VERITAS experiment. 
A somewhat surprising result from our study is that the information of the 
various 7-hadron separation methods is largely uncorrelated and combining 
the information from several parameters, one obtains a much more powerful 
hadron suppression. We conclude with a summary and a discussion in Sect. 6. 



2 Simulation Details 



We used the Grinnell-ISU-Utah (GrISUU) air shower and detector simulation 
package^ that combines the KASCADE air shower simulation code [24] with 
the calculation of the Cherenkov light emitted by the air shower and the 
simulation of the detector response. We generated 480,000 vertically incident 
7-ray-initiated air showers over the energy range from 30 GeV to 10 TeV, 
distributed in energy according to a differential power law spectrum dN/dE oc 
E~^ with photon index F = 2.5. The 7-rays were assumed to originate from 
a point source located at the center of the field of view of the telescopes and 
were simulated over a circular area of 350 m radius. We simulated 1,930,000 
proton-initiated air showers over the energy range from 100 GeV to 20 TeV 
with a power law index of 2.7 and with arrival directions uniformly distributed 
over a 4° radius circular area centered on the field of view of the telescopes. 
The proton-initiated air showers were generated over the same area as the 7- 
ray-initiated ones. All simulations assumed that the experiment is located at 
an altitude of 1.8 km above sea level. The GrISUU code traces the incoming 



http://www.physics.utah.edu/gammaray/GrISUl 



Cherenkov photons through the mirror and hght cone geometry and simulates 
the response of the photomultiphers (PMTs) and the digitization of the signals 
with the VERITAS 500 MHz flash analog to digital converters (FADCs). It 
simulates single pixel triggers, and the trigger of a telescope. We set the pixel 
trigger threshold to five photoelectrons. A telescope triggers if three pixels 
fire within 10 nsec and satisfy the pattern trigger requirement of any three 
adjacent pixels. The telescope array triggers if three telescopes trigger within 
50 nsec. The night-sky noise was set to 4.2x10^^ photoelectrons Hz m~^ sr~^. 

While the condition of three telescopes triggering in coincidence results in a 
higher energy threshold than a two-telescope trigger condition, it may opti- 
mize the sensitivity of the array for the majority of sources detectable with 
VERITAS. For sources with a soft GeV spectrum, a two-telescope trigger 
may give a higher sensitivity. When comparing the performance benchmarks 
of different analysis methods (e.g. angular and energy resolution, 7-hadron 
separation capability), one should be careful to specify the trigger condition. 
The methods discussed below can certainly also be applied to data taken with 
a hardware two-telescope coincidence condition, imposing a software require- 
ment of a detected image in at least three telescopes. The choice between a 
two-telescope and a three-telescope trigger condition depends also on other 
considerations, as e.g., the dead time of the data acquisition system. It may 
be the best choice to decide on the trigger condition after taking some real 
data with various trigger conditions, rather than to rely entirely on simula- 
tions. In the following, we will sometimes use events where at least three or 
at least four telescopes triggered. All figures shown in this paper assume the 
three-telescope trigger condition. 

The events were analyzed with the "eventdisplay" package developed at the 
University of Leeds. The PMT charge is determined by fitting a sample trace 
to the PMT traces and numerically integrating the best fit function. The code 
uses pedestal events to determine the average pedestal variance (RMS of the 
pixel charge in the absence of an air shower related signal). The pedestal vari- 
ances are determined by the rate of night sky background photons and the 
level of electronic noise. The analysis characterizes PMT charges in units of 
the pedestal variance of each pixel. The image cleaning process consists of set- 
ting all amplitudes to zero for which neither the "image" (the signal amplitude 
exceeds four times the pedestal variance) nor the "border" (a pixel adjacent 
to an image pixel with a signal amplitude exceeding two times the pedestal 
variance) condition is met. 

In the following we present distributions of several air shower parameters. 
Note that all parameter distributions discussed in this section have been de- 
rived from the true parameter values used in the Monte Carlo simulations. 
Distributions of reconstructed parameters will be discussed in the following 
sections. If not stated otherwise, we show distributions for all events that trig- 



gered the experiment and that had at least three telescopes with an image 
(more than two pixels surpassing the image threshold). 

The left panel of Fig. 1 shows the distribution of air shower cores of the 
triggered 7-rays (the points where the air shower axes intersect the telescope 
plane). Air showers basically trigger when their cores fall within the air shower 
array or less than ~150 m away from the outer three telescopes. The right 
panel of Fig. 1 gives the distribution of the distances of air shower cores from 
the center of the telescope array for both triggered 7-rays and protons. The 
core distributions smoothly approach zero at 350 m from the central telescope, 
showing that we simulated air showers over a sufficiently large area to produce 
a realistic air shower dataset. In Fig. 2 we present the distribution of the en- 
ergies of triggered events for both 7-rays and protons. The simulations cover 
sufficiently broad energy ranges, as the distributions approach zero close to 
the boundaries of the simulated range. The differential detection rate per log- 
arithmic energy interval peaks at about 100 GeV for photons and 550 GeV for 
protons. We see that proton-initiated events have a substantially higher energy 
threshold than 7-ray-initiated events. The latter is a well known consequence 
of (a) hadronic showers channeling a fraction of their energy into muons and 
neutrinos via vr+'^-production, and (b) hadronic showers being more irregular 
than purely electromagnetic showers and thus being more unlikely to generate 
multi-pixel and multi-telescope coincidences. 

Fig. 3 presents the effective detection area for 7-rays. It rises proportional to 
E"^'^ in the threshold region (<100 GeV) and proportional to E^''^'^ at higher 
energies (>300 GeV). 

The distance of the arrival directions of triggered proton events from the center 
of the field of view is shown in Fig. 4. At ~4° from the center of the field of 
view, hardly any protons trigger the telescopes. 



3 Reconstruction of Shower Axis 



The reconstruction of the air shower direction and core is based on the simple 
fact that the major axes of the images are projections of the shower axis. 
In the field of view of the telescopes the major axes intersect at the point 
corresponding to the arrival direction of the primary particle. In the reference 
frame of the telescopes, all ellipse-like images point away from the shower 
core (the intersection of the shower axis and the telescope plane) and the core 
can again be found at the intersection point of the major axes [11,12]. If the 
cosmic 7-ray source is a point source at a well-known position in the sky, the 
air shower core can be determined with very high statistical accuracy. One 
can use the fact that the lines pointing from the location of the source to the 



centroids of the air shower images point to the air shower core [13]. These hnes 
can be reconstructed with a higher statistical accuracy than the major axes 
of the elhpse images. If not stated otherwise, we will assume in the following 
that the source is a point source of known location. 

In practice, finding the best reconstruction algorithm means combining the 
information from all telescopes in an optimal way to minimize the statisti- 
cal error of the direction and core estimates. Here we treat the problem as 
a x^-minimization problem. First the direction is determined by finding the 
point that minimizes the weighted squared distances to the major axes. Sub- 
sequently, shower core is reconstructed by minimizing the weighted squared 
distances to the lines that go through the location of the source and the cen- 
troids of the images. Compared to averaging over intersection points from pairs 
of telescopes, this approach has the advantage that the information about the 
location and orientation of all telescopes enters the estimate at the same time 
rather than sequentially. As usual in a x^"fit, the weights used for combining 
the information of all images should be inversely proportional to the "squared 
statistical error" associated with each image. We experimented with three dif- 
ferent weighting schemes. The simplest ones use constant weights or weights 
inversely proportional to the parameter size, the sum of the counts of the cor- 
responding image. Furthermore, we scrutinized the mean distances between 
the major axes (or lines between the source direction and image centroids) and 
the true direction and core as function of several parameters, i.e., size, width, 
length, width/length, and distance of the telescope to the shower axis r. The 
mean distances indeed scale approximately proportional to 1/size. After scal- 
ing the errors with 1/size, we still find a rather strong residual dependence on 
width (not on width/length as one may naively expect as the width/length- 
values characterize the ellipticities of the images). We determined an empirical 
weighting function that depends on the size, width and length parameters. 
Fig. 5 presents the angular and core resolutions achieved with three weighting 
schemes. While the more sophisticated weighting schemes improve on the an- 
gular and core resolutions, the improvement is very small. Weighting according 
to the third scheme, we obtain an angular and core resolution (63 % value) 
of 0.22 ° and 7.5 m, respectively. The lack of a substantial improvement may 
stem from the fact that the images in the same shower tend to have the same 
quality, minimizing the effect of the weights. The results agree well with those 
of an earlier study performed for the HEGRA Cherenkov telescope array [13]. 

The performance improves if high-quality events are selected. Using only 
events with a reconstructed shower core within 150 m from the central tele- 
scope and with a reconstructed energy exceeding 300 GeV we obtain angular 
and core resolutions of 0.1° and 4.1 m, respectively. Figure 6 shows the angu- 
lar resolution for different software energy thresholds. The resolution improves 
with increasing energy threshold and above 1 TeV we obtain a resolution of 
0.05°. The cut on the core location is critical for obtaining the good results. 



Without it, more and more events far away from the telescope array trigger the 
array. At the hardware trigger threshold the cut excludes 24% of the photon 
initiated events from the analysis. 

For extended sources, we do not know the arrival direction of the photons a 
priori, and the major axes have to be used to determine the shower core loca- 
tion. Rather independent of the applied weighting scheme, the core resolution 
deteriorates from 7.5 m (4.1 m) for the point source analysis method to 37 m 
(17 m) for the extended source analysis method (the values in brackets have 
been determined for a software energy threshold of 300 GeV). 



4 Energy Reconstruction 



Together, the atmosphere and the Cherenkov telescopes constitute a fully 
active calorimeter with sparse sampling. The energy of the primary particle 
is roughly proportional to the Cherenkov light intensity measured with the 
telescopes. In the following we discuss only the reconstruction of the primary 
energy for 7-rays. We implemented a simple energy reconstruction method 
(see also [19,25,26]). We determined the median and 90%- width- values of the 
logarithm of the size parameters as function of the primary 7-ray energy E and 
distance r from the shower axis. For each telescope with a telescope trigger, 
an energy estimate is determined by inverting the lookup table. The energy of 
the primary particle is determined by averaging the energy estimates from all 
telescope with a telescope trigger. We weight the estimates from each single 
telescope proportional to one over the square of the statistical uncertainty 
on the estimate. Combining in this way the energy estimates from different 
telescopes gives an energy resolution CT\■c,(E,^c/Eu^^^) of 0.28 for the analysis of 
point sources and 0.4 for the analysis of extended sources. We focus in the 
following on the analysis of point sources. 

In order to improve on the energy estimate we take into account that the 
Cherenkov light intensity depends on the height of the shower maximum ^ . 
Given the stereoscopic air shower data, the height of the shower maximum 
can be determined. We use the ratio between the distance of a telescope to 
the shower axis r and the angular distance dist between the centroid of the 
image and the reconstructed shower direction as an estimate of the shower 
maximum. For a constant shower maximum we expect that r and dist are 
roughly proportional to each other and thus that r/dist= const. The closer 
the shower maximum is to the observation level, the further away the image is 
from the center of the field of view and the larger is the dist value for a given 



^ A shower reaches the "shower maximum" when the number of electrons/positrons 
that emit Cherenkov light reaches a maximum. 



distance from the shower axis, and the smaller is the r/distTaiio. The left panel 
of Fig. 7 shows the logarithm of the ratio between the reconstructed energy and 
the true energy as function oir/dist. Indeed, we see that the simple algorithm 
severely underestimates (overestimates) the energy for showers with a shower 
maximum high (low) in the atmosphere. We add a correction that depends on 
the r/dist ratio and the size measured in a telescope. The dependence on the 
size-parameter takes into account that the correction depends somewhat on 
the primary energy of the inducing particle. The right panel of Fig. 7 compares 
the energy resolution achieved with and without the correction. The correction 
improves the energy resolution cri^{Ercc/Etruc) of the air shower array from 0.28 
to 0.22. The performance improves if high-quality events are selected. Using 
only events with a reconstructed shower core within 150 m from the central 
telescope and selecting events with a reconstructed energy exceeding 300 GeV 
we obtain an energy resolution o'i^{Ercc/Etruo) of 0.15. 



5 Gamma-Hadron Separation 



One of the main strengths of Cherenkov telescope experiments is the large 
collection area on the order of 10^ m^. As a consequence, they achieve un- 
equaled sensitivity for 7-ray observations on short time scales. However, for 
longer integration times their sensitivity is limited by fluctuations in the rate 
of background events and thus increases only with the square root of the ob- 
servation time. At energies of 100 GeV the air shower background consists of 
Cosmic Ray hadrons and Cosmic Ray electrons. We only discuss the suppres- 
sion of the first. Based on today's technology, electron initiated air showers 
can only be suppressed based on their isotropic arrival direction. In the follow- 
ing we study only proton-initiated air showers. Protons are expected to make 
up about 75% of the Cosmic-Ray initiated background. The other 25% are 
mainly produced by Cosmic Ray He-nuclei. We expect that proton-initiated 
air showers resemble 7-ray-initiated air showers more closely than air showers 
initiated by heavier nuclei. 



Normalized Width 



We discuss five methods to distinguish between 7-ray and hadron-initiated 
air showers (see also the excellent discussion in [27,28]). The first approach 
analyzes the width of the air shower images perpendicular to the major axes. 
Hadron-initiated air showers show significantly "wider" images than 7-ray- 
initiated showers owing to the transverse momentum inherent in hadronic 
interactions. The transverse momentum originates in the non-negligible ki- 
netic energy of the nucleons inside hadronic nuclei and the quarks inside the 
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nucleons. Similar to Aharonian et al. (1997) [11], we use the width parameters 
measured in all telescopes to derive a 7-hadron separation parameter. Based 
on the Monte Carlo dataset, we derive a lookup table of expected median 
width values, Wm, and the 90%- widths of the distributions, Wqq, as functions 
of the size of the image and the distance r of the telescope from the shower 
axis. Here the median and 90%-width are used rather than the average and 
the RMS to reduce the impact of outliers in the width distributions. Here and 
further below, we have used half the dataset for optimizing the method, and 
the other half to measure its performance. 

Given the w^{r,size) and WQo{r, size)-va\ues from the lookup table, a "nor- 
malized width" value is computed by: 



w = — — X 

-' 'trig 



1: '""""'""""•'•""■> 



Wgo[ri,size 



where the sum runs over all iVtrig telescopes with a telescope trigger, widths 
and sizci are the width and size values of the image found in the i**^ telescope 
with a trigger, and r[ is the distance of the i*'^ telescope from the shower axis. 
The distribution of the normalized width parameter for photons and protons 
is shown in the left panel of Fig. 8. Accepting events with a normalized width 
value below a certain cut value Wcnt for the analysis, will accept a fraction e^ 
of 7-ray-initiated events and a fraction of ep of proton-initiated events. The 
Q-factor of a 7-hadron separation method is defined as 



Q = ^ (2) 



where e^ and Cp are the fractions of photon and hadron initiated events surviv- 
ing the cut. The Q-factor resembles the improvement in signal to noise ratio 
achieved with a cut, assuming that the noise is dominated by background 
fluctuations. The values e^, ep and the Q-factor are given as function of the 
cut value widthcut in the right panel of Fig. 8. We obtain an optimal Q-factor 
of 1.5 for Wcut = 0.3. Using the more restrictive condition that four rather 
than three telescopes triggered, improves the performance of the normalized 
width cut to a maximum Q-value of 2.38 for Wcut = -0.1. The performance of 
the normalized width cut, as well as that of the 7-hadron separation methods 
described in the following paragraphs is summarized in Table 1. 
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Agreement Between Different Telescopes Regarding the Shower Direction and 
Core 



Krennrich & Lamb (1995) [29] and Hillas (1996) [27] suggested to use the 
deviation of the major axes from the reconstructed shower direction and/or 
core location as means to distinguish between photon and hadron-initiated 
events. The wider and more irregular hadron-initiated showers tend to produce 
images that do not all point to the arrival direction and core location. We use 
chi-square values Xdir ^^'^ Xcore ^f the direction and core fits, respectively, 
to differentiate between photons and hadrons. We normalize the chi-square 
values according to the number of telescopes participating in the fit. We find 
that the optimal cuts in the chi-square values Xdir ^-nd Xcorc result in hadron 
suppressions with Q-factors of 1.34 and 1.69 (three-telescope trigger condition) 
and 1.53 and 2.3 (four-telescope trigger condition), respectively. The cut in 
Xcore is ^^^y powerful, and achieves a comparable 7/hadron separation as the 
cut in the normalized width parameter. 



Lateral Cherenkov Light Distribution 



The fourth method takes into account that photon and hadron initiated 
air showers exhibit a different lateral Cherenkov light distribution, and that 
Cherenkov light pool is more homogeneous for photon initiated air showers 
than for hadron initiated air showers. We make use of these differences by 
comparing the energy estimates derived from the information of individual 
telescopes with the energy derived from the information of all telescopes. 
Hereby, all energy estimates assume that the primary particle is a photon. 
We define a reduced chi-square value as follows: 



iVtrig i i^i ^hi(Ei) 



The sum runs over all telescopes with a trigger, E^a is the energy estimate 
derived from all telescopes, Ei is the energy estimate derived from the i*'' 
telescope, and o"in(£;.) is the estimated error on ln(_Ei). Optimizing the cut- 
value, we find that the x^-cvX gives maximum Q-factors of 1.30 and 1.60 for 
the three-telescope and four-telescope trigger conditions, respectively. 
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Making Use of the Temporal Information 



Another consequence of the more regular nature of the electromagnetic show- 
ers compared to hadronic showers, is that their Cherenkov light front is nar- 
rower in time. The VERITAS array of Cherenkov telescopes is equipped with 
fast FADCs to read out the pixel signals. As fifth 7-hadron separation method, 
we compute for each pixel the time at which the signal rose to 50% of its peak, 
and compute a reduced chi-square value Xtimc ^^ ^^e root mean square (RMS) 
value of these rise times divided by the degrees of freedom. For the three- 
telescope and four-telescope trigger conditions, a cut in Xtimc gives Q-factors 
of 1.09 and 1.22, respectively. 



Combining the Information of Various Methods 



Scatter plots between two of the five 7-hadron parameters discussed above 
show that the parameters are not strongly correlated with each other. Cut- 
ting in several parameters should result in an improved performance. We apply 
the likelihood ratio formalism to combine the information from different ap- 
proaches. We use the 7-ray and proton distributions of the cut parameters as 
probability density functions (PDFs). The 7-hadron separation parameter is 
then the logarithm of the ratio of the probability that the event is a photon 
divided by the probability that it is a proton. We define the parameter: 

N 

Ai = E [i^(^") - HPH] (4) 

i=l 



where the sum runs over the A^ 7-hadron separation methods to be used, 
and P-^ and P^" are the probabilities that photon and proton-initiated air 
showers produce the observed value of the i*^ parameter, respectively. While 
this approach should give the optimal results if the parameters are indeed 
completely uncorrelated, we also compute a second parameter which might 
work better if there is some correlation between the parameters: 

A2 = min (ln(i^^) - In(i^P)) (5) 

1=1. .N 



Basically, we identify an event as a proton-initiated event if one of the used 
parameters strongly favors the proton over the photon hypothesis. The per- 
formance of several combinations of the separation parameters are given in 
Table 1. For the three-telescope coincidence condition, we find that combining 
the information from all five 7-hadron separation methods improves on the 
best Q-factor 1.7 of a single method, to a Q-factor of about 2.6 (see Fig. 9). In 
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the case of the four-telescope trigger condition, the parameter w alone gives 
Q =2.4. Combining the information from w and Xcorc improves the perfor- 
mance to a Q-factor of 3.6. In all cases, there is little difference between using 
Equ. (4) or Equ. (5) for combining the information. 

Please note that we get very comparable results and almost identical 7-hadron 
separation Q-factors for the analysis of point sources and extended sources. 

Finally, we would like to mention that imposing the four-telescopes trigger 
condition reduces the number of detected 7-ray and cosmic ray events. Taking 
into account the numbers of detected events, the achievable angular resolu- 
tions, and the best Q-factors, we compute that the signal to noise ratio of a 
point source detection is by a factor of 1.9 better for the events taken with 
the four-telescope trigger condition than for all the events taken with the 
three-telescope trigger condition. 



6 Summary and Discussion 



We have used a dataset of simulated photon and proton air showers to study 
the performance of simple event reconstruction methods. Furthermore, we 
have explored different schemes to improve on the performance of the meth- 
ods. We have shown that the direction and core reconstruction are insensitive 
to details of the weighting of the major axis of the individual telescopes. The 
energy resolution benefits from correcting for the height of the shower max- 
imum. For the analysis of point sources, we get angular, core and energy 
resolutions of 0.22° (0.1°), 7.5 m (4.1 m), and 22% (15%), respectively. The 
first set of numbers applies to all events that produce a three-telescope trigger. 
The set of numbers in brackets applies to events with a shower core within 
150 m from the central telescope and with a reconstructed energy exceeding 
300 GeV. In general, our performance estimates for the VERITAS experiment 
agree well with the experimentally verified performance of the H.E.S.S. exper- 
iment [30]. 

We have compared different 7-hadron separation methods and have shown 
that the information from several parameters can be combined based on the 
likelihood ratio approach. Compared to a cut in a single 7-hadron separation 
parameter, a cut on the likelihood ratio improves the Q-factor from 1.7 to 2.6 
(three-telescope trigger) and from 2.4 to 3.6 (four-telescope trigger). 

It is instructive to compare the relative importance of the hadron and electron 
backgrounds before and after applying the 7-hadron separation cuts. Above 
30 GeV the BESS energy spectrum for protons impinging on the atmosphere 
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is [31]: 



— ^ = 9.6 X 10-9(^/1000 GeV)~2.7^j^-2g-ig^-iQgy-i ^g^ 

dE 



The 1 GeV-100 GeV electron spectrum measured by HEAT is [32]: 

dF 

^ = 1.2xl0-^x (^/lGeV)-^x 
dE 

(l + {E/5 GeV)2-3) "^ cm-2 s"^ sr^^ GeV^^ (7) 



In Fig. 10 we show the proton and electron trigger rates as function of the true 
and reconstructed shower energy before and after applying a 7-hadron sepa- 
ration cut. We used here the three-telescope trigger condition and a 7-hadron 
separation cut in Ai(w, Xdir^ Xcore; xD- -^^ proton-initiated showers produce 
less Cherenkov light than purely electromagnetic showers, the reconstructed 
proton energies are on average lower than the true ones. For electrons, there is 
not such a systematic shift. Notwithstanding this dramatic effect, the hadron 
background still dominates over the electron-background - even after applying 
the gamma-hadron separation cut. The result shows that excellent 7-hadron 
suppression continues to be of utmost importance for VERITAS and other 
similar experiments. 
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Table 1 

Q-factors, photons and proton acceptance and cut value for the 7-hadron separation 

methods described in the text. 
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Fig. 1. The left panel shows the distribution of the core locations of the trig- 
gered 7-ray-initiated air showers in the plane of the telescopes. For reference, 
the circles give the location of the 4 VERITAS telescopes. The right panel 
shows the distribution of the distances of the air shower cores from the cen- 
tral telescope Tl for 7-ray (solid line) and proton (dashed line) initiated air 
showers that triggered the experiment. 
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Fig. 2. Distribution of the primary energy of photons (solid line) and protons 
(dashed line) that triggered the telescope array. The photon distribution peaks 
at ~100 GeV and the proton distribution at ~550 GeV. 
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Fig. 3. Effective area for detecting cosmic 7-rays as function of primary en- 
ergy. Above tlie tliresliold region (>100 GeV), tlie effective area increases only 
slowly with energy. 
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Fig. 4. Distribution of the angular distance between the arrival directions of 
all protons that triggered the experiment and the center of the field of view 
of the cameras. 
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Fig. 5. The left panel shows the point spread function for all photon-initiated events 
that triggered the telescope array. The different line styles show different weighting 
algorithms. The solid line shows the results for constant weights; the dashed line 
shows weights proportional to the size parameter, and the dotted line shows a weight 
that depends on the size, width and length parameters. All three weighting schemes 
give very similar results. For the best method (dotted line), we get an angular 
resolution (63% value) of 0.22°. The right panel shows the distributions of the error 
in the core location for all 7-ray-initiated events. The different line styles correspond 
to the same weighting schemes as describes for the left side. For the best method 
(dotted line), we get a core resolution (63% value) of 7.5 m. 
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Fig. 6. Angular resolution (63% value) as function of threshold energy. The 
reconstructed energy was used to impose the threshold energy cut. Only events 
with reconstructed core locations within 150 m from the central telescope were 
used here. The angular resolution substantially improves above the hardware 
energy threshold. 
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Fig. 7. The left panel shows the logarithm of the ratio between reconstructed 
and true air shower energy as function of the ratio between the distance from 
a telescope and the dist parameter measured in that telescope. The latter ratio 
depends somewhat on the height of the shower maximum. The higher (lower) 
in the atmosphere the shower develops, the larger (smaller) the ratio. One can 
see that the reconstructed energy is overestimated (underestimated) for show- 
ers developing low (high) in the atmosphere. The right side compares the en- 
ergy resolution achieved with correction (dashed line) and without correction 
(solid line) for the height of the shower maximum. After (before) correction we 
get an energy resolution of 22% (28%). The plots include all events for which 
the reconstructed energy exceeds 100 GeV and for which the reconstructed 
shower core less than 200 m away from the central telescope. 
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Fig. 8. The left panel shows the distribution of the normalized width pa- 
rameter for photon (solid line) and proton (dashed line) initiated events 
(three-telescope trigger condition). The right panel shows the photon accep- 
tance e^ (dashed line) and proton acceptance ep (dashed-dotted line) as func- 



tion of a cut on the normalized width parameter 
shows the Q-factor. Here it peaks at Qmax = 1-5. 
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Fig. 9. Same as Fig. 8 for the cut \i{w, Xdn^ Xcorcy xD- The maximum Q-factor 
is 2.60. 
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Fig. 10. The top three curves show the rate of proton events found in a circular 
bin of 0.25° radius at the center of the field of view of the telescopes. The 
top three curves show the proton rate (i) before application of the 7-hadron 
separation cut as function of the true proton energy (solid line), (ii) before 
cut as function of the reconstructed energy (dashed line), and (iii) after cut as 
function of the reconstructed energy (dotted line). The lower three lines show 
the same for the electrons rather than protons. One can see that the proton 
background dominates over the electron background, even after applying the 
7-hadron separation cut. We used here the cut \i{w, Xdir? Xcoro xD- 
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